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We study the entanglement dynamics of a system consisting of a large number of coupled harmonic oscillators 
in various configurations and for different types of nearest neighbour interactions. For a one-dimensional chain 
we provide compact analytical solutions and approximations to the dynamical evolution of the entanglement 
between spatially separated oscillators. Key properties such as the speed of entanglement propagation, the 
maximum amount of transferred entanglement and the efficiency for the entanglement transfer are computed. 
For harmonic oscillators coupled by springs, corresponding to a phonon model, we observe a non-monotonic 
transfer efficiency in the initially prepared amount of entanglement, i.e., an intermediate amount of initial 
entanglement is transferred with the highest efficiency. In contrast, within the framework of the rotating wave 
approximation (as appropriate e.g. in quantum optical settings) one finds a monotonic behaviour. We also study 
geometrical configurations that are analogous to quantum optical devices (such as beamsplitters and interferom- 
eters) and observe characteristic differences when initially thermal or squeezed states are entering these devices. 
We show that these devices may be switched on and off by changing the properties of an individual oscillator. 
They may therefore be used as building blocks of large fixed and pre-fabricated but programmable structures in 
which quantum information is manipulated through propagation. We discuss briefly possible experimental re- 
alisations of systems of interacting harmonic oscillators in which these effects may be confirmed experimentally. 
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I. INTRODUCTION 

Quantum information processing requires as a basic ingre- 
dient the ability to transfer quantum information between spa- 
tially separated quantum bits, either to implement a joint uni- 
tary transformation or, as a special case, to swap quantum 
information between the qubits. For a transfer over larger 
distances it is usually imagined that some stationary qubits, 
for example in the form of trapped ions inside an optical 
resonator, are coupled to a quantized mode of the electro- 
magnetic field that propagates between the spatially separated 
cavities |ll0|3l01- Other specific realisations are possible, 
but the basic principle always relies on the use of some con- 
tinuous degree of freedom between the qubits and their manip- 
ulation by external fields. While this appears to be the most 
realistic mode of transport over long distances, one may con- 
ceive other modes over shorter distances. Instead of using a 
quantum field one may study the possibilities offered by a dis- 
crete set of interacting quantum systems. This might involve 
spin degrees of freedom or infinite dimensional 

systems such as harmonic oscillators JlollTlll . In the present 
paper we explore the dynamics of entanglement in a chain of 
coupled harmonic oscillators llT2lfl3lfl4ll . Apart from its obvi- 
ous relevance to quantum optical systems including photonic 
crystals, such a model also describes phonons in a crystal and 
we therefore hope that the results presented here will also have 
applications in condensed matter systems as well. 

This paper is organised as follows. In Section 2 we present 
the basic physical models and their Hamiltonians. We use 
analytical tools from the theory of Gaussian states in continu- 
ous variable systems where some rapid development has been 
achieved recently (for a tutorial overview see, e.g., Ref. HI). 
We briefly reiterate those results that will be employed in the 



present investigation. The entanglement properties of a sys- 
tem of harmonic oscillators in the static regime have been 
studied in some detail fl2il . Section 3 will then present the 
basic equations of motion in compact form for two types of in- 
teractions, namely (a) harmonic oscillators coupled by springs 
and, resulting from this, (b) a model which corresponds to 
a rotating wave approximation as is appropriate in a quan- 
tum optical setting. Section 4 employs these equations for the 
propagation of entanglement along a chain of harmonic oscil- 
lators which might be realized by coupled nano-mechanical 
oscillators JTsIIIct or optical cavities. The time evolution of 
the entanglement between a pair of oscillators is given ana- 
lytically in a compact form. Properties such as the speed of 
propagation, the amount of entanglement and the transfer ef- 
ficiency are then obtained from these expressions. In Section 
5 we present results utilizing the equations from Sections 3 
and 4. Firstly, we study a method for the creation of entangle- 
ment in such a system that does not require detailed control of 
the interaction strength between individual oscillators but only 
the ability for changing the interaction strength globally lllll . 
The influence of imperfections such as finite temperatures or 
randomly varying coupling constants on such a scheme are 
studied. We also consider the propagation of some initially 
prepared entangled state along the chain. Surprisingly, for 
harmonic oscillators coupled by springs we observe a non- 
monotonic transfer efficiency in the initially prepared amount 
of entanglement, i.e., an intermediate amount of entanglement 
is transferred with the highest efficiency. Conversely, in the 
rotating wave approximation, the transfer efficiency is mono- 
tonic. While most of these results assume a position indepen- 
dent and stationary coupling we also show that with carefully 
chosen position dependent coupling the transfer efficiency in 
this system may be increased to unity. Finally we study geo- 
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metrical configurations that are analogous to quantum optical 
devices such as beamsplitters and interferometers and observe 
characteristic differences when initially thermal or squeezed 
states are entering these devices. We show that these devices 
may be switched on and off by changing the properties of 
an individual oscillator and may therefore be building blocks 
of large fixed but programmable structures. In Section 6 we 
summarize the results of this paper and suggest possible ex- 
perimental realisations of systems of harmonic oscillators in 
which these effects may be confirmed. 



11. MODELS AND METHODS 

In this section we present the systems under considera- 
tion, namely coupled harmonic oscillators, together with the 
Hamiltonians that describe the various models for their inter- 
action. We will restrict our attention to Hamiltonians that are 
quadratic in position and momentum operators. This will be 
crucial for the following analysis as it permits us to draw on 
the results and techniques from the theory of Gaussian contin- 
uous variable entanglement. The most important results from 
this theory will be reviewed here briefly. 



A. The physical models 

The general setup consists of a chain of AI coupled har- 
monic oscillators, where the coupling is assumed to be such 
that the corresponding Hamiltonian is at most quadratic in po- 
sition and momentum. We will number the harmonic oscil- 
lators from 1 to M with periodic boundary conditions such 
that the {M + l)-th oscillator is identified with the 1-st. The 
choice of periodic boundary conditions yields exact and com- 
pact analytical solutions since we can employ normal coordi- 
nates straightforwardly. A similar approach is less successful 
in the non-periodic boundary case. In addition, we allow for 
the existence of a distinguished decoupled oscillator with in- 
dex which will be a convenient notation for some of the later 
studies. Arranging the position and momentum operators in 
the form of a vector 

R = (90, 91, ■ • ■,qM,Po,Pi, ■ ■ ■ ,Pm) (1) 
we can then write the general Hamiltonian in the form 



energy is p1/2. As a consequence, both the potential 
matrix and the kinetic matrix are diagonal and identi- 
cal, namely V ^ T = 1m+i, where 1La/+i denotes the 
M + 1 by A/ + 1 identity matrix. 

(b) (Oscillators coupled by springs) If neighbouring oscil- 
lators (except for the 0-th oscillator) are coupled via 
springs that obey Hooke's law, the Hamiltonian is given 
by 
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(3) 



where c denotes the coupling strength and we have used 
periodic boundary conditions i.e., (7j\/-|_i ~ qi. Keeping 
in mind that we wish to leave the oscillator with index 
uncoupled, the potential matrix becomes 



V = 



1 
l + 2c 
-c 1 







-c 



I 

c 

2c -c 

c 1 + 2c 








-c 








■■■ -c 

-c l + 2c -c 
-c l + 2c 



(4) 

while the kinetic matrix is given by the identity matrix 
T = 1La/+i- 

(c) (Oscillators in rotating wave approximation) An inter- 
action that provides simpler dynamics is obtained via 
the rotating wave approximation in quantum optical 
systems. Indeed, if we define annihilation and creation 
operators 



a = -^{q + ip), 



(5) 



then we observe that the interaction in case (b) includes 



terms of the form a 



i.e., interaction terms for 



which both harmonic oscillators are being excited si- 
multaneously. Such terms are not energy conserving, 
and in quantum optics they are usually be neglected 
in the framework of the rotating wave approximation 
(RWA). Following this practice amounts to considering 
the following Hamiltonian 



H 



1 



V 
T 



1 



M 



= ^^qiVi^qj +PiTi]Pj, (2) 

•y = l 



where V is the potential matrix and T the kinetic matrix. We 
will consider three basic settings for which we now provide 
the matrices T and V explicitly. In all these cases we assume 
that the oscillators in the chain are all identical with a mass 
TO = 1 and eigenfrequency lo ~ 1. 

(a) (Uncoupled oscillators) If the oscillators are not cou- 
pled to each other, then the potential energy of the k- 
th oscillator is simply given by qf./l while its kinetic 



RWA 
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\- CLk+ial 



(6) 



In terms of position and momentum operators this can 
be written as 



9o +Po 



1 



Qk +Pk 



fe=l 



+ I (9fc+i - Qkf + ^ {Pk+i - Pkf , (7) 
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so that in this case both the potential and the kinetic 
matrix are given by 



T = V 
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(8) 



Note that the matrix V can be conceived as the adjacency 
matrix of a weighted graph G(u, e) encoding the interaction 
pattern between the systems in the canonical coordinates cor- 
responding to position. Vertices of the graph correspond to 
physical systems, i.e., the individual harmonic oscillators, 
whereas the weight associated with each of the edges quan- 
tifies the coupling strength flTll . The main diagonal corre- 
sponds to loops of the weighted graph. This intuition is in 
immediate analogy to graph states for sp in syste ms with an 
Ising interaction between the constituents 1 18, 19,^20] and can 
be useful in the study of more complex geometries. Before we 
study the dynamical properties of these systems, we provide 
in the following subsection a brief overview over the main 
technical tools that we are going to employ. 



B. Analytical tools 

Analysing the entanglement properties of infinite dimen- 
sional systems such as harmonic oscillators is generally tech- 
nically involved unless one restricts attention to specific types 
of states. Indeed, in recent years a detailed theory of so-called 
Gaussian entangled states has been developed. As we will em- 
ploy some of its basic results in the subsequent analysis we are 
providing a brief review of some useful results. A more de- 
tailed tutorial review can be found, e.g., in Ref. fl(ill . and more 
technical details concerning Gaussian states can be found in 
Ref. (2H. 

The relevant variables in the analysis of harmonic oscilla- 
tors are the canonical operators for position and momentum. 
Let us assume a system with n harmonic oscillators. As stated 
above it is convenient to arrange these in the form of a vector 
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,Pn) ■ 



The characteristic feature distinguishing the quantum har- 
monic oscillator from its classical counterpart is the canon- 
ical commutation relation (CCR) between position and mo- 
mentum. Employing the vector R these can be written in the 
particularly convenient form [Rj, Rk] = i<Tj,k where the real 
skew-symmetric block diagonal 2n x 2n-matrix a, the sym- 
plectic matrix, given by 



In 

-^n 



(9) 



assuming units where = 1, and fc = 1, a choice that will 
be adopted for the rest of this paper Instead of referring to 
states, i.e., density operators, one may equivalently refer to 
functions that are defined on phase space. While there are 
many equivalent choices for phase space distributions, for the 
purposes of this work it is most convenient to introduce the 
(Wigner-) characteristic function. Using the Weyl operator 
= e*^ '^^ for ^ G R^", we define the characteristic func- 
tion as 



Xp(6=tr[pM/d- 



(10) 



The state and its characteristic function are related to each 
other according to a Fourier- Weyl relation. 



P - 
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^ / d^"exp(-Ow-«. (11) 



Gaussian states are exactly those states for which the charac- 
teristic function Xp is a Gaussian function in phase space 
That is, if the characteristic function is of the form 



(12) 



As is well known, Gaussians are completely specified by their 
first and second moments, d and 7 respectively. As the first 
moments can be always made zero utilizing appropriate local 
displacements in phase space, they are not relevant in the con- 
text of questions related to squeezing and entanglement and 
will be ignored in the following. The second moments can be 
collected in the real symmetric 2n x 2?! covariance matrix 7 
defined as 



7j,fc 



2Retr[p(4- - {Rj)p){Rk - {Rk)p)] ■ (13) 



With this convention, the covariance matrix of the n-mode 
vacuum is 7 = l2„. 

As the covariance matrix encodes the complete informa- 
tion about the entanglement properties of the system, we will 
use it in order to quantify the amount of entanglement be- 
tween two groups of oscillators. There is no unique way to 
quantify entanglement for mixed states, and several differ- 
ent measures grasp entanglement in terms of different oper- 
ational interpretations. For the purposes of this work we set- 
tle for the logarithmic negativity which is comparatively easy 
to compute and possesses an interpretation as a cost function 
iHllllllllllil. Given two parties, A and B, the logarith- 
mic negativity is defined as 



N{p)^log,\\p' 



(14) 



where p-^^ is the state that is obtained from p via a partial 
transposition with respect to system B and 1 1 . 1 1 1 denotes the 
trace-norm. As we focus attention on Gaussian states which 
we characterize via the covariance matrix 7 rather than the 
density matrix p, we need to provide a prescription for the 
evaluation of the logarithmic negativity directly in terms of the 
covariance matrix. To this end, it is important to note that on 
the level of covariance matrices the transposition is reflected 
by time reversal which is a transformation that leaves the po- 
sitions invariant but reverses all momenta, q ^ q, p —p. 
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The partial transposition is then the application of this time- 
reversal transformation to a subsystem, i.e., one party. Let us 
now consider a system made up of m + n oscillators, where m 
oscillators are held by party A and n oscillators by party B. 
Applying time reversal to the n oscillators held by party B, 
the covariance matrix will be transformed to a real symmetric 
matrix 7^^ given by 



where 



P = l, 



!,„ 
-1L„ 



(15) 



(16) 



The 2n x 2n-matrix 7^^ is the matrix collecting the second 
moments of the partial transpose of p. The logarithmic 
negativity is then given by 



N = 



- log2(min(l, 



l7,l)), 



(17) 



where the 7^ are the symplectic eigenvalues of 7^^ . For a gen- 
eral covariance matrix, jj arises in the diagonalization of 7 
using symplectic transformations, i.e., transformations S that 
preserve the CCR so that SaS'^ = a. The resulting diagonal 
matrix is the Williamson normal form of a covariance matrix 
whose diagonal elements are the symplectic eigenvalues. It 
is sometimes useful to know that the symplectic eigenvalues 
can be obtained directly without explicit diagonalization of 
the matrix as the positive square roots of the usual eigenval- 
ues of — (T7CT7 |24]. 

For all Hamiltonians that are quadratic in the canonical po- 
sition and momentum operators the ground state is an impor- 
tant example of a Gaussian state. For a Hamilton operator of 
the form 



V 
T 



R 



(18) 



we find that the covariance matrix of the ground state is given 
by 



7 



which reduces to 



(19) 



(20) 



when T = 1L„. If on the other hand, as for the interaction 
in case (c), we have T = V, then the ground state is given 
by 7 = 1„ 1„, which is the same as the ground-state of n 
non-interacting harmonic oscillators. 

The primary aim of this work is the investigation of the dy- 
namical properties of the system of harmonic oscillators and 
the evolution of entanglement properties under such dynam- 
ics. The dynamics of the covariance matrix under a Hamil- 
tonian quadratic in position and momentum operators can be 
obtained straightforwardly from the Heisenberg equation 



dt 



X{t) = i[H,X]. 



(21) 



For our time-independent Hamiltonian Eq. ( I18> . this leads to 
the covariance matrix at time t as 



lxx{t) "fx pit) 
1px(t) 7pp(t) 



exp 



X cxp 



T 
~V 

-V 
T 



Ixx IXP 
IPX IPP 



(22) 



Equipped with these tools we can now proceed to the analysis 
of the entanglement dynamics of the harmonic chain. 



III. THE EQUATIONS OF MOTION 

In two separate subsections, we provide the explicit solu- 
tions to the equations of motion for the two coupling models 
(b) and (c) that will be investigated both analytically and nu- 
merically in the remainder of the paper 



A. Harmonic oscillators coupled by springs 



This model is characterized by a Hamiltonian of the form 

M 

2 



1 9 

^Spring ^ T^Y^k+Pl+C^ - Qkf 



(23) 



fc=i 



Note that for the moment we neglect the decoupled additional 
0-th oscillator In the following we will provide an explicit 
form for the equations of motion for the canonical positions 
and momenta in the Heisenberg picture. To this end we can 
diagonalise this Hamiltonian by introducing the normal coor- 
dinates 



1 \ ^ 27Tin 

— > e 

M 

m— 1 

1 



(24) 



This leads to 



iJspring = \Y, (PsPl + C^fQ.Ql) (25) 

where cj^ = 1 + 4csin^(7rs/Af). Here we have used the fact 
that Qu = Q'Lu^Pu = P-u- introduce the annihilation 
operators 



-.(ujsQs+iPl), 



so that the Hamiltonian takes the form 

M 



Spring 



alas + 



s=l 



(26) 



(27) 
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In the Heisenberg picture the annihilation and creation oper- where we have defined the useful functions 
ators evolve according to as{t) = e^"^=*a5(0) and al{t) = 
g+jwst^t (^0). Separating the real and imaginary parts, we get 



Qs{t) = Q.(0)cos(w3t) + — Pj(0)sin(cj,t), (28) 

Us 

Ps{t) = -lOsQI{0)sui{loJ) + Ps{0)cos{ujst) . 

Substituting these into Eq. J24> . we obtain the time evolution 
for the original position and momentum operators 

Qnit) = ^((?.(0)/,_„(t)+p,(0)g,_„(<)), (29) 

r=l 
M 

Pn{t) = ^(q,(0)/,_„(i)+Pr(0)/.-„(t)), 



r=l 



J 



1 " 

fk{t) = J^YI cos(a;„i), 



(30) 



m—l 

M 



9k( 



In the entire paper, dots will denote time derivatives. Defin- 
ing the covariance matrix elements to be (once we ignore the 
displacements (gi)p) 



79„<j„ = 2Re tr[/9g„9„ 
we find their values at time i as 



(31) 



X] (/r-nW/s-m(i)79r«s + 5r-n W^s-m (t)7p.p« + /r-« (O^s-m (079.P, + 5r-n (OA-™ (07p.g, ) : (32) 
A/ 

{jr-n{t)h-m{i)lq^q, + 5r-n (t)/s-m (i)7p.p, + fr-n{t)fs-rn(t)-/q^p^ + Qr -n{t) j s-m{t)lp,q}j , 

',s=l 
M 

^ ^ (^fr — n(t)fs — m(j')^qrqs fr — n(i)fs — m(t)^PrPs ~^ fr — n(t)fs—m(j')^qj.ps ^ fr — nii) fs — m(j')^prqs^ ^ 



where the 7 on the right hand side are the initial values of the 
covariance matrix elements. 



B. Hamiltonian in the rotating wave approximation 

This model is characterized by the Hamiltonian 

Hrwa = 2 51 ^fe + 2 ^ ^''-'^ ^ 2 ^ ^''''^ ' 
fe=i 

(33) 

and in the following, we can carry out the analysis along the 
same lines as in the previous subsection. Again, employing 
the normal coordinates Eqs. (I24> we obtain 

M 

HnwA = 2 E {P^Ps + QsQl) (34) 



s=l 



where we now have 



n,^l + 2csin^ 



M 



(35) 



Introducing the annihilation and creation operators 



±{Qs+^P^^), at ^_L(Qt„,p^) (35) 



the Hamiltonian assumes a form 



M 



(37) 



In the Heisenberg picture the annihilation and creation oper- 
ator then evolve in time as as{t) ~ e~'^^s*as(0) and dj(t) — 
e+*f^s*Q,t(o). Separating again real and imaginary parts we 
obtain 



Qs{t) = Q.(0)cos(f|2t) + pt(o)sin(0^t), (38) 



Transforming back to the original position and momentum op- 
erators we find 



M 

<in{t) = ^(g,(0)F,_„(t)+p,(0)G,_„(t)), (39) 

M 

Pn{t) = ^(-(Zr(0)G.-„(i)+p.(0)F,_„(t)), 
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where we have defined the functions and Gk as 

A/ 1 

^'^■(^) = 77 E cos(17,2„t), G,.(t) = TT E -^"^ ^^(^'^t) ■ (40) 



r/i— 1 m— 1 



Note that these functions are slightly simpler than the corresponding ones in Eqs. (|30} as they lack the frequency ftg in the 
denominator. The covariance matrix elements vary in time as 
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— Gr-n{i)Gs-77i{'^)^pr-qs^ 
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r,s— 1 


^Gr-n{t)Gs-ni{t)^qrqs ~^ 




Gr—n(i)Fs — m{t^'~fqrPs ~ 


Fr—n(j'^Gs — 'ni(j')'^Pj~qs^ ; 



(41) 



IV. PROPAGATION OF ENTANGLEMENT ALONG THE CHAIN 



We would like to investigate the capacity of the harmonic chain for transmission of quantum information. The clearest 
signature for the ability to transmit quantum information and coherence is the proof of the ability to transmit one constituent 
part of an entangled pair of oscillators through the chain. To analyze the situation we require the equations of motion for the 
covariance matrix. Let us assume the existence of a distinguished oscillator which is entirely decoupled from the others. We 
imagine that at time t = this oscillator is prepared in a two-mode squeezed state with the first oscillator of the chain while all 
other oscillators are assumed to be in their respective ground state (assuming no interaction). The initial conditions then read 

Iqoqo = Iqiqi = IpoPo = IpiPi = cosh(r), = -7poPi = sinh(r), (42) 

7g.9. = IpsPs = 1 for all s > 1, 
79.P. . 



The 0-th oscillator for the interaction via springs will obey a free time evolution which is given by (using Eqs. ( I28> and noting 
that go = Qo and j5o = Po) 

go(0 = go(0) cos(tJo^) + , po(0 = -'^0<7o(0)sin(cJo0 +Po(0)cos(a;ot)- (43) 

Similarly, for the RWA interaction 

9o(i) = qo(0)cos(f]2t)+po(o)sin(r!2t), po(i) = -<7o(0) sin(172t) + ^^(o) cos(^]2t) . (44) 

Note, however, that they correspond to local unitary rotations on the 0-th oscillator which do not affect the entanglement between 
this oscillator and the remaining ones. To simplify the expressions we will therefore omit this time evolution in the following. 
Again we will treat the two types of interactions described by iJspring and Hrwa separately. The next two subsections establish 
the analyutical expressions for the time evolution of the entanglement between the 0-th oscillator and the n-th oscillator 

A. Harmonic oscillators coupled by springs 

In this case, employing the special form of the initial conditions for the system Eqs. J42> . the elements of the covariance 
matrix describing the 0-th and the ?i-th oscillator at time t are given by 

79090 (0 = 7popo(0 = cosh(r), 79o9n(0 = sinh(r)/„_i(i), 
7poP,. (t) = ~ sinh(r)/„_i (t), jq^p^ (t) = 0, 
JqoP^it) = sinh(r)/„_i(t), Jq„po{t) = - sinh(r)g„_i(i), 

(45) 
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and 

M 



s=l 

M 



Iq^pAt) = (cosh(r) - l)(/„_i(t)/„_i(t) +g„_i(t)/„_i(t)) +3„_.(t)/„_,(t), 

s=l 

M 

Ip^pAt) = (cosh(r)-l)(/2_,(t) + /2_i(t)) + ^(/2_,W + /,t,W). 



s=l 

In the limit of a chain of infinite length, i.e., when M — > oo, we find 

Iqoqoit) = Ipopoit) cosh(r), 79o9„(0 = 'IpoP^i*) = sinh(r) /„_ i (t) , (47) 
7gopo(0 0. IqoP^it) = sinh(r)/„_i(i), 

7goP„(*) = sinh(r)/„_i(t), Jq„po{t) = - sinh(r).g„_i (i), 

and 

7<?„g„(i) = (cosh(r)-l)(/2_i(t)+.g2_^(i))+a„„(t)+rf„,^(i)^ (48) 
lq,^P^{^) = (cosh(r) - l)(/„_i(t)/„_i(t) +g„-i(t)/„-i(0) + &,m(0 + e„„(t), 
7p„p„(i) = (cosh(r) - + fl_^{t)) + c„„(t) + a„„(t). 



Here we have employed the definitions C ~ c/(l + 2c) and Q, = \/l + 2c, and introduced the functions 

1 /"^ / \ I r ^ N sii^ (mVl - 2Ccos( 



fk{t) = — [ dcf) cos{k(j)) cos (vityjl — 2Ccos^) , gA:(i) = — / dcj) cos{k<f)) 



OVl - 2Ccos( 



and 











COS ^ 




7^ Jo 





a„„(t) = - 1 + -/ d^cos 2rJtVl-2Ccos(/) , (49) 



bnnit) 



— [ d(j) n^/l- 2Ccos0sin (217^1 - 2Ccos(t>) , 



2 o2 



Cnn{t) — — -Z— — 2CC0S(/)) COS f2ili-\/l — 2(^ COS ^) , 

2 27r Jo ^ ^ 

, , , 1 1 1 r cos (20 VI - 2CCOS0) 

«rm(0 = -^===-— / ^ 



At) 



2 VI + 4c 27r7o 02(l-2Ccos( 
1 sin (20V1 - 2C cos 



27r Jo OVl - 2Ccos( 



While the above set of equations determines the time evolution exactly, they do not provide very much direct insight into the 
dynamics of the system. In the following we will show however, that we can obtain very good and compact approximations 
to the above exact solution in terms of elementary functions. While the following derivation is not rigorous in the sense that it 
does not provide error bounds, a numerical comparison between the approximations and the exact results shows the impressive 
quality of the approximate solution. As a first step, we expand the functions fk and gk to first order in 



I r If 

fk{t) ~ — d(j) cos {k(j)) cos {ilt {I — ( cos (f>)) , gk{t) — — / rf(/)Cos(A;(/)) sin (Ot(l — ^cos^)) , (50) 
Jo 71- Jo 

and also drop a factor of 1 — 2Ccos(/) in gk- In the following we will employ Bessel functions which satisfy the relations 

oo 

cos(a;coss) ~ Jojx) + 2 J2n{x) coa(2ns) cos(n7r), 

n=l 

oo 

sin(a;coss) = 2 J2„+i(x) cos((2n + l)s) cos(n7r) . (51) 
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On using trigonometrical addition theorems one finds that in first order 



9k{t) ^ JkiC^^t) sm{nt-^ 



(52) 



A further crucial approximation replaces the time-dependent quantities a„„(t), bnn{t), c„„(t), dnn{t), e„„(i) by their time 
averages, i.e.. 



ann{t) i 



lim — 

T-,00 T 



lim — 



lim — 



dtann{t) 



dtCnn{t) 



1 

Imi — / dtbnn{t) = 0, 

r-»oo 1 Jq 



(53) 



0= 



1 dnn{i^ 



lim — 



dtd„„{t) 



dtCnnit) = 0. 



With all these approximations we finally obtain 

79090 W - 7poPo(<) - cosh(r), 

7909„ (t) 
Iqopo (^) 
7goP„ (i) 

79„g„(i) 

7p„p„(i) = (cosh(r)-l)J;;_i(Cm) + l + c 



(54) 



7pop„ (^) = sinh(r) J„_i(CrJt) cos {vit 

7<7„pn it) 

(cosh(r)-l)j2_i(Cl)t) 



7r(n — 1) 



sinh(r) J„_i sin 

1 1 1 



7r(n — 1) 



2 2ym^' 



Numerical comparisons in later sections will show that these approximate solutions are very good approximations when r is not 
too small. We have so far collected all the elements of the covariance matrix involving the 0-th and the n-th oscillator. Since we 
will investigate the entanglement properties between the two oscillators, we can trace out the rest of the oscillators which leaves 
us with the covariance matrix of the reduced system comprising only two oscillators. Employing the ordering (qo,j5oj QnjPn) 
we find 



7Spr 



A D 
B 



where 



A 



B 



D 



cosh(r) 
cosh(r) 



(cosh(r) - l)Js{tf + 1/2 + (1 + Ac)-^'^/2 

(cosh(r) -l)Js{tf 

sinh(r) Js(i) cos $s(i) - sinh(r) Js(t) sin$s(t) 
- sinh(r) Js(i) sin<i>s(t) - sinh(r) Js(0 cos $,s(i) 



(55) 



(56) 



1 +c 



where we have used the abbreviations ^s{t) ^ fit — ^{n — l)/2 and Js{t) = Jn-i{C^t). From this explicit form for the 
covariance matrix of the 0-th and n-th oscillator we can now determine the symplectic eigenvalues as solutions of the polynomial 



T]^ - (dct(yl) + dct(B) -2dct{D))if + det(7spring) = 0. 



(57) 



The solution is given by 



^Spring 



(2/1 - y2) , 



(58) 
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where we have 



and 



2 + c + (1 + c)(l + 4c)-' - (5 + 2c + (1 + 4c)-5) Jf + iJ^ 
+Jl{i + 2c + (1 + 4c)"5 _ 4j2) cosh(r) + (1 + J^f cosh(2r), 



2/1 



2/2 = \ w + yf 



(59) 



= -8[2jf + (! + (! + 4c)-3 _ 2J2) cosh(r)][j2 + (1 + c - J^) cosh(r)]. 



(60) 



Note that we have dropped the time argument in Js {t) to make the expressions appear more compact. The logarithmic negativity, 
finally, is in this approximation given by 



^Spring(0 = -log2(niin(|7ysp,.i„g|,l)). 
Note that the other symplectic eigenvalue is greater than one. 



(61) 



B. Hamiltonian in rotating wave approximation 

For the Hamiltonian i?RWA one proceeds along very similar lines, i.e., taking the limit A/ ^ oo to find 

1 



Fk{t) = 
Gkit) = 



1 



)COs(fc0) COs[f]RwA*(l - CrWA COS0)], 

)cos(fc0) sin[12|\vAi(l - Crwacos0)], 



(62) 



with r^RWA = vT+c and Crwa = c/(l + c), which are exactly 

Fn{t) = Jn{ct) cos{nly^j^t ~ 7m/2), G„(i) = Jn{ct) sm{nly^j^t ~ Tm/2) . (63) 
we find the covariance matrix elements to be 

Iqoqoit) = 7popo(0 = cosh(r), (64) 

7go«„W = -7poP„W = sinh(r)K-i(i), 
7goP„W = 7g„poW -sinh(r)G„_i(i), 
7,„«„ W = 7p„p„ W - (cosh(r) - l){F^,_,{t) + Gl_,{t)) + 1. 

Note that now the terms corresponding to Unnit) , bn7iit) , Cnn{t) , dnnit), and enn{t) do not need to be approximated as all 
time-dependencies cancel each other off conveniently in the expressions as opposed to the spring case. Again we can write the 
covariance matrix of the reduced system comprising only the 0-th and the n-th oscillator employing the ordering {qa,pa, Qn ,Pn), 
and find 



7RWA 



A D 
B 



with 



A 



B 



D 



cosh(r) 
cosh(r) 

(cosh(r) - 1) JRWA(t)' + 1 

(cosh(r)-l)jRWAW' + l 

sinh(r) Jrwa(0 cos $rwa(0 - sinh(r) Jrwa(0 sin $rwa(0 
- sinh(r) Jrwa(<) sin $rvv^a(0 - sinh(r) jRWA(i) cos $RWA(i) 



(65) 



(66) 



where we have used the abbreviations $rwa(^) = ^rwa^ ^ sions are very similar to those for the first interaction type. 
n{n — l)/2 and JRWA(i) — Jn-i{ct). Note that these expres- 
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The symplectic eigenvalues are given by the solutions of 

(det(^)+dct(B)-2det(L»))r?^WA+det(7RWA) = 0. 

(67) 

This gives rise to 

'^RWA = ^ (-1 - ^2) , (68) 

with 

zx = (1 + Jrwa) cosh^(r) + 2 J^^a sinh^(r) 

+2^wa(1 - ■^rwa) cosh(r) + (1 - J^sNpf . 
Z2 = (69) 

where 

V = -4 ((1 - J^wa) cosh(r) + J2^a)' ■ (70) 

Again we have dropped the explicit time dependence in Jrwa 
for brevity of notation. The logarithmic negativity is then 
given by 

A^RWA(i) = - log2(min(|?7RWA|, 1)) ■ (71) 

Having prepared all the analytical work we need, we can now 
proceed to investigate the entanglement dynamics of the har- 
monic chain. 



V. STUDY OF THE ENTANGLEMENT DYNAMICS OF 
THE HARMONIC CHAIN 

This section studies numerically and analytically, various 
aspects of the entanglement dynamics of the harmonic chain 
with various initial states and geometrical arrangements for 
both types of interactions. We begin by briefly revisiting the 
effect of the spontaneous creation of entanglement which is 
obtained when the interaction strength between the oscilla- 
tors is globally changed suddenly LI 1.1 . This effect can only 
be observed in a model in which the oscillators are coupled 
via springs. In the rotating wave approximation this effect 
does not occur as both the coupled and uncoupled chains have 
an identical ground state. Then we move on to consider the 
propagation of entanglement through a harmonic chain. Sur- 
prisingly, for harmonic oscillators coupled by springs corre- 
sponding to a phonon model, we observe a non-monotonic 
transfer efficiency in the initially prepared amount of entan- 
glement, i.e., an intermediate amount of entanglement is trans- 
ferred with the highest efficiency. In the rotating wave approx- 
imation, the transfer efficiency is monotonic though equally 
surprising. We provide analytical expressions for the propa- 
gation speed of the entanglement through the chain together 
with approximate analytical expressions for the transfer effi- 
ciency. We also study the influence of imperfections such as 
finite temperatures and varying coupling constants. Finally we 
study different geometrical configurations that are analogous 
to quantum optical devices such as beamsplitters and interfer- 
ometers and observe characteristic differences when initially 



thermal or squeezed states are entering these devices. We pro- 
pose ways in which these devices may be in and out of action 
thereby allowing for the creation of pre-fabricated quantum 
networks that can be programmed by external switches. 

A. Spontaneous creation of entanglement 

We consider all the harmonic oscillators to be in the ground 
state and initially uncoupled, i.e.. 

We suddenly switch on the interaction at time t = {) and ob- 
serve the dynamical evolution of entanglement. We do not ob- 
serve any entanglement with the i?RWA interaction (because 
we have not included a mechanism to produce a spontaneous 
excitation, that is, we do not have terms such as a\a\^i) so 
we will exclusively deal with iJspring- We are limiting our 
scope to numerical results because the analytical expressions, 
while they can be provided, are too complicated to yield any 
insight. Note also that the approximations used below in the 
squeezed state case cannot be applied here. A typical exam- 
ple of the time evolution of entanglement is shown in Figure 
n In an open chain of length 30 we study the time evolution 
of the entanglement between the first and the last oscillator 
We observe no entanglement for a finite time until at a time 

0.1 2 I : i i i 1 




50 100 150 200 250 

t 



FIG. 1 : The time evolution of entanglement between the first and the 
last oscillator in an open chain of length 30 when all the oscillators 
are initially in the ground state. The coupling strength c = 0.1 has 
been chosen for the solid line (c = 0.15 for the dashed line). 

one first encounters a build-up of entangled between the two 
oscillators. This time is approximately given by 



ie the time is approximately linear in n, the separation of 
the oscillators. It should be noted that <o is half as large as 
the time that is required for a perturbation at the first oscil- 
lator to travel to the n-i\\ oscillator This suggests that the 
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origin of the entanglement between the l'** and the n-th os- 
cillator arises from the interaction of those oscillators exactly 
half-way in between. Their entanglement is generated by the 
initial sudden switching on of the interaction and then prop- 
agates through the chain. This idea will be further corrob- 
orated in the later subsections when the propagation of pre- 
prepared entanglement is considered. Furthermore one finds 
that the dependence of the maximal degree of entanglement as 
measured by the logarithmic negativity is approximately pro- 
portional to n~^^^ for large n until reaches values of about 
10~^ when it begins to drop quite rapidly to vanish entirely. 
It should be noted that for any parameters of the model there 
exists a finite n such that the state of the first and the n-th 
oscillator is separable for all times. For the coupling value 
of c = 0.1 this will happen for n = 20000. Therefore, the 
value of largest separation is very large indeed for reasonable 
coupling strengths. To appreciate that this is a somewhat sur- 
prising behaviour, it should be contrasted with the entangle- 
ment structure in the ground state of the chain. Then it can 
be shown that for any chosen coupling strength two distin- 
guished oscillators are never in an entangled state, unless they 
are immediately neighbouring |,1 J|. 

Here we have studied the special case of the entanglement 
between the endpoint of an open chain. It should be noted 
that this is a particularly favourable configuration. For a given 
distance between oscillators one always obtains the largest 
amount of entanglement when one places them at the oppo- 
site ends of an open chain. Two oscillators in a very long chain 
with the same distance and with positions well away from the 
ends of the chain will lead to considerably smaller amounts of 
entanglement. Indeed, the amount of entanglement will differ 
by approximately a factor of 4. This discrepancy is due to the 
fact that at the ends of the chain the oscillators possess fewer 
neighbours with which they can become entangled. As we dis- 
card all oscillators other than two any entanglement with other 
oscillators will deteriorate the entanglement between the dis- 
tinguished oscillators. While this does not explain the factor 
of 4 quantitatively, it gives an intuitive picture for the decrease 
of entanglement that will discussed in more detail later on. 



interaction. The thermal equilibrium state is given by 

2 



2 



(74) 



The equations of motion are then 



2 


pu/T _ 


1 


2 






1 


2 






1 



(t) + dnrnit)), (75) 



7p„P,„ (i) = (l + pc^/T _ I ) (Crmi(^) + a„,„(t)). 

Appropriate values for the temperatures have to be obtained 
from experiment in the particular context under consideration, 
but it appears possible nowadays to achieve ratios T/cj <C 1 
in different physical systems (note that we have taken ft = 1 
and fc = 1) such as nano-mechanical oscillators. We consider 
again an open chain instead of a closed ring. This renders the 
analytical treatment more difficult but makes no difference for 
the numerics. The entanglement depends little on the tem- 
perature as long as the mean thermal photon number is well 
below unity. Figure |2l shows the temperature dependence of 
entanglement evolution. We observe that down to tempera- 




Time 



B. Thermal state case 

In the previous section we have studied the entanglement 
dynamics in an environment that is at zero temperature. This 
is reflected by the fact that the initial state of our harmonic 
oscillators is assumed to be the vacuum. We will now move 
one step further towards a more realistic description by set- 
ting the initial state as a thermal equilibrium state. As with 
the ground state case, we do not establish entanglement for 
the i/RWA interaction as a thermal state can be represented 
as a mixture of displaced vacuum states which do not lead to 
any spontaneous entanglement in the rotating wave approxi- 
mation. Therefore we shall again only consider the i/spring 



FIG. 2: The time evolution of entanglement between the first and 
the last oscillator (with open boundary conditions). We have fixed 
c — 0.1, the chain consists of 30 oscillators, and x — cu/T for 
X > 10 (solid line) and x = 6 (dashed line). 

tures corresponding to x = 10, the entanglement evolution 
is almost exactly the same as in the ground state case. Only 
when X < 10 do we see an effect of a finite temperature. Even 
for x = 6, we still have a significant portion of entanglement 
albeit a small delay in the arrival of entanglement. A more 
realistic scenario is dealt with in Ref. iHHi . 

Decoherence mechanisms may be included without leav- 
ing the harmonic setting. Often, the high temperature limit 
of Ohmic quantum Brownian motion is appropriate with an 
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independent heat bath for each oscillator in the limit of negli- 
gible friction and under the assumption of product initial con- 
ditions. Such a decoherence mechanism can be accounted for 
by adding terms of the form iIztIEsIi 



<[qn, [qn,p]] 



(76) 



to the idealized unperturbed generator of the dynamical map 
for each of the oscillators, where the real number ^ specifies 
the decoherence time scale. However, in cases where product 
initial conditions are inappropriate or unrealistic, decoherence 
may still be modeled using, for example, time-convolutionless 
projection operator techniques ll29il . In small systems, non- 
product conditions may be incorporated by explicitly append- 
ing heat baths to each of the oscillators with a linear coupling 
13011 . according to Hamiltonians of the form 



H„ 



(77) 



(i) 

with real numbers ^j, where the denote the canonical co- 
ordinates corresponding to position of the i-th oscillator of 
the j-th heat bath consisting of m oscillators. Assuming a 
particular form of the spectral density, the coupling strength 
to the finite heat baths may be chosen in a way that is consis- 
tent with empirically known values for energy dissipation. Of- 
ten, Q-factors are approximately known for resonators, which 
quantify the number of radians of oscillations necessary for 
the energy to decrease by a factor of 1 /e. Hence, on the basis 
of these Q-factors, the appropriate coupling may be evalu- 
ated. Figure |3] shows the influence of decoherence in case of 
an open chain with the same parameters as in Figure^for an 
Ohmic heat bath, i.e., = jA/m, where A > is the cut-off 
frequency of the environment modes. One finds that the cre- 
ated entanglement by suddenly switching on the interaction is 
surprisingly robust against decoherence within this model. 




Time 



FIG. 3: The same as Figure Q without Ohmic dissipation and de- 
coherence (solid line), and with Ohmic dissipation corresponding to 
Q = 10000 (dashed line) and Q = 1000 (dotted line). 



C. Entanglement transport through the harmonic chain 

In the previous section we considered the case where no en- 
tanglement was present in the initial state of the system. En- 
tanglement emerged as a consequence of a sudden change in 
the coupling constant between neighbouring harmonic oscil- 
lators. In this subsection we are going to investigate a differ- 
ent situation and consider the transmission of entanglement 
through a one dimensional chain. To this end we initialise 
two harmonic oscillators in a two mode squeezed state. We 
assume that one of these oscillators is decoupled from the rest 
and give it the index 0, while the other oscillator, with index 
1, forms part of a chain of harmonic oscillators with near- 
est neighbour interaction. The remainder of the chain starts 
out in the ground state corresponding to zero temperature (as- 
suming no interactions). By evolving the initial state, we ex- 
pect the entanglement to travel along the chain so that with 
increasing time more and more distant oscillators will be en- 
tangled with the 0-th harmonic oscillator. There are a number 
of free parameters that can be varied: the coupling strength 
c, the amount of initial entanglement quantified by the two- 
mode squeezing parameter r, the time t and the position of 
the oscillator to be entangled with the 0-th oscillator. In order 
to simplify the analytical work, we shall be dealing with the 
limit M ^ oo. 

While in the previous section the interaction in the rotating 
wave approximation does not lead to the spontaneous creation 
of entanglement, here it allows for the propagation of entan- 
glement. We give an example for the time evolution of the 
logarithmic negativity between the 0-th and the 30-th oscil- 
lator for both interactions in Figure |3 For both interactions 
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FIG. 4: The time evolution of entanglement between the 0-th and the 
30-th oscillator in a chain with 80 oscillators with periodic boundary 
conditions. We have chosen c = 0.1 and r = 0.8 in both cases. The 
entanglement propagates slightly faster in the rotating wave interac- 
tion. 



we obtain qualitatively the same behaviour but we observe 
that under the RWA interaction the entanglement propagates 
somewhat faster but as expected this difference decreases with 
decreasing coupling constant c. Another difference is the fact 
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that the entanglement under the RWA interaction does not ex- 
hibit the small-amplitude oscillations that the interaction due 
to harmonic oscillators coupled by springs exhibits due to the 
existence of counter-rotating terms of the form Sfc 0.^+1. The 
propagation of the quantum entanglement can be seen even 
more clearly in Figure |5l where for a ring composed of 40 
oscillators and a coupling constant of c = 0.1 the time evo- 
lution of the logarithmic entanglement between an uncoupled 
oscillator and the n-th oscillator is shown when initially the 
uncoupled oscillator and the 1-st oscillator are coupled. One 
observes that with increasing time more and more distant os- 
cillators are becoming entangled. Entanglement propagates 
both clockwise and anti-clockwise around the ring. After a 
sufficiently long time it becomes important that the ring has a 
finite size and the two counter-rotating 'entanglement waves' 
meet at the opposite end of the ring and we observe some en- 
tanglement enhancement. Both Figures |3 and [S] suggest that 




FIG. 5: For a ring of 40 oscillators and a coupling constant of 
c = 0.1 the time evolution of the logarithmic entanglement between 
an uncoupled oscillator and the n-th oscillator is shown when ini- 
tially the uncoupled oscillator and the 1-st oscillator are in a two- 
mode squeezed state with two-mode squeezing parameter r = 0.8. 
With increasing time more and more distant oscillators are becom- 
ing entangled. Entanglement propagates both clockwise and anti- 
clockwise around the ring. After a sufficiently long time the two 
counter-propagating 'entanglement waves' meet at the opposite end 
of the ring and we observe some entanglement enhancement. 



entanglement can be distributed to distant oscillators. It will 
therefore be interesting to study the efficiency for this transfer 
when we vary the amount of entanglement provided initially 
by varying r. In particular, we will be interested in the first lo- 
cal maximum in the amount of entanglement Nf as quantified 
by the logarithmic negativity. We separate the study of the en- 
tanglement transfer efficiency for the two interactions as they 
exhibit distinctly different behaviours. We begin with the in- 
teraction describing oscillators interacting via springs. Figure 
|6l shows the amount of entanglement at the first local maxi- 
mum as squeezing is varied. We observe the remarkable fact 
that for large initial entanglement, the value of N j saturates. 
We can obtain an analytic expression for the saturation value 



0.15 




01 23456789 10 

Squeezing parameter r 

FIG. 6: For oscillators interacting by springs the graph shows the 
amount of entanglement at the first local maximum as squeezing is 
varied. We have fixed c = 0.1 and n = 30. 



as a function of c and n by taking the r 00 limit. Thus, we 
approximate ?7spring — (Syi) of Eq. (I58> and, discarding 
all the terms except cosh^(r) ^ e^''/4 and cosh(2r) = e^''/2, 
we obtain 

2 ^ (-2 JLx + (1 + 4c)-i/^ + 1) {-Jl,^ + c + 1) 

where J^ax = 0.6748851(n - ly^/^ ^ is the value of 
the first maximum of the 71-th Bessel function of the first 
kind. This substitution provides a very good approximation 
as the first local maximum of the logarithmic negativity co- 
incides with that of the n-th Bessel function. Therefore we 
find iVsat = -log2(|??Spring|) which is shown in Figure H 
We observe that the saturation value decreases for both in- 




FIG. 7: A graph to show how the saturation value A*'sat for the loga- 
rithmic negativity varies with the coupling strength c and the position 
n of the second oscillator. 



creased coupling strength and increased distance. While the 
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latter is intuitive the former might be surprising as one could 
have thought that it is advantageous to increase the coupling 
strength to facilitate the transfer of entanglement but this in- 
tuition is clearly contradicted by Figure Even more strik- 
ingly the entanglement vanishes entirely when the interaction 
strength becomes too high. We believe that this is due to the 
fact that the initial entanglement disperses across several os- 
cillators and will discuss this at the end of this subsection. 

If we translate these findings into an entanglement transfer 
efficiency defined by 



cff 



(79) 



then we observe that this efficiency exhibits a non-monotonic 
behaviour Indeed, we observe a maximum in the efficiency 
as shown in Figure|S]for the same parameters as in Figure|5] 
We have not yet found a convincing and intuitive explanation 
for the occurrence of a maximum in the transfer efficiency 
js^l- In fact we will shortly see that the phenomenon of a 
non-monotonic transfer efficiency is absent in the RWA in- 
teraction. The surprising implication of this non-monotonic 
behaviour of the transmission efficiency is that it is advanta- 
geous to transmit entanglement in intermediate size portions 
rather than in one very large packet. 




Squeezing parameter r 



FIG. 8: The efficiency of transmission as squeezing is varied for c — 
0.1 and n = 30. Surprisingly, the transmission efficiency exhibits a 
maximum at a finite initial entanglement. 



Let us now consider the same problem of the entanglement 
transfer efficiency in the RWA interaction. Indeed, we can 
show that there is still saturation in the amount of entangle- 
ment that can be transmitted and we find the value of the sat- 
uration to be (after taking r oo) 



Nsat = - l0g2 



J- 

1 + j; 



2 

max 



(80) 



where again J„iax = 0.6748851(n - 1)-^/^ |131|] is the value 
of the first maximum of the n-th Bessel function of the first 
kind. Note that this expression is independent of the cou- 
pling strength c. Unlike the case of oscillators interacting with 



springs there is no maximum in the efficiency for the RWA in- 
teraction as can be seen clearly in Figure |5] Indeed, for large 
r the efficiency is tending to zero while for r approaching 
the efficiency tends to Toff = Jmax- 

Since the efficiency is not equal to unity for both interac- 
tions, the question arises as to where the rest of the entangle- 
ment is located? The most obvious place is to search in the 
neighbourhood of the n-th oscillator. Since we always deter- 
mine the entanglement between individual oscillators we ig- 
nore many others that have interacted with it and have thereby 
become entangled. Any entanglement between these two os- 
cillators will therefore deteriorate as they are being entangled 
with other oscillators that we choose to ignore. This viewpoint 
is corroborated by determining the entanglement between the 
0-th oscillator and a whole group of neighbouring oscillators 
instead of a single one. The result of this can be seen in Fig- 
ure [TO] This shows the change in the amount of entangle- 
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FIG. 9: The efficiency of the RWA interaction versus the two-mode 
squeezing parameter r,n — 30 and c — 0.1. 



ment as we increases the number of oscillators in the second 
group. The graph supports the idea that the missing entan- 
glement between the 0-th and the n-th oscillator is due to the 
creation of entanglement between the n-th oscillator and its 
neighbours because we start to recover entanglement as we 
compute the entanglement between the 0-th oscillator and the 
neighbourhood of oscillators surrounding the ?i-th oscillator. 
This spread of entanglement is not dissimilar to a dispersion of 
the energy of a wavepacket as it experiences different group 
velocities. However, the effect on the entanglement can be 
considerably stronger as the energy will only decrease linearly 
with the width of the wave-packet while the entanglement can 
drop much more rapidly and become zero at finite spreading. 



D. Speed of entanglement propagation 

We have found that the propagation of half the two mode 
squeezed state through the chain takes a finite time as the 
entanglement between the 0-th and the n-th oscillator is ex- 
actly zero for a finite time interval (see for example Figure|3i. 
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unlike the case of interaction via springs, the efficiency un- 
der the rotating wave approximation does not decrease as we 
increase c because its efficiency is independent of c. 




FIG. 10: The change in the amount of entanglement as the number 
of oscillators of system B centered around the 20-th oscillator is in- 
creased from 712 = 1 to 712 = 21. We have a ring of 70 oscillators 
with squeezing parameter r — 0.8. Increasing the number of neigh- 
bours increases the amount of entanglement available. This supports 
the idea that the loss of entanglement is due to dispersion. 



After a certain amount of time, the two oscillators in ques- 
tion become entangled and the logarithmic negativity reaches 
a temporary maximum. We are able to determine this time 
analytically for both types of interactions that we are consid- 
ering. To make the analysis tractable, we consider an infinitely 
long chain. We find that the first maximum of the logarithmic 
negativity coincides with the first maximum of a Bessel func- 
tion Jm{x). We know the position of this maximum occurs 
at .T = m + 0.8086165771^/3 HI. Noting that m = n - 1, 
X ~ (^ilt = ct/ \/l + 2c for the Hooke's law interaction and 
a; = ct in the RWA interaction, we obtain 



n- l + 0.8086165(n- 1)1/3 

prmg - ^ : 

n- 1 + 0.8086165(71- 1)1/3 



^RWA 



(81) 



We observe that the time that is required for entanglement to 
be established between the 0-th and the ri-th oscillator is a 
function of the coupling strength and the position ri. For large 
separations n, it quickly becomes linear in 7i and as expected, 
the larger the coupling c the faster entanglement is established. 
We also see that the RWA interaction produces faster entan- 
glement since tspring/iRWA = Vl + 2c. As 71 — 1 is the 
separation of the 1-st and the 7i-th oscillator we can define the 
speed of propagation to be 



'''^P™*^ VTT2^ (1 + 0.8086165(71 - l)-2/3) VTT2^' 



E. Optimization of entanglement transfer and generation 

In the previous section we have studied the entanglement 
transfer along a chain of identical harmonic oscillators as this 
will be the situation that is most easily implemented exper- 
imentally. However, we observe that the transmission effi- 
ciency decreases with distance. One might expect that one 
can improve this efficiency by tuning the couplings and the 
eigenfrequencies of the harmonic oscillators suitably. Indeed, 
in this section we will show what can be achieved in this more 
general setting. For simplicity we consider the task of trans- 
mitting one half of a two-mode squeezed state from one end 
of an open chain to the other 

We assume as usual one decoupled harmonic oscillator with 
index and a chain of length M through which the other half 
of the two-mode squeezed state is transmitted. Perfect trans- 
mission from one end of the chain to the other is possible 
in the rotating-wave interaction of nearest neighbours if we 
choose the interaction strength 



Vn,n+1 = Vn+l,n = C\/ n{M - Tl) 



and 



1, 



(83) 



(84) 



with the real number c being sufficiently small in order for V 
to be positive. The choice of the diagonal elements being all 
equal to 1 is equivalent to the requirement that we choose the 
eigen-frequencies w„ of the uncoupled oscillators as 



l-cV(Af-l), 



(85) 



n 



1 - cy/n{M - n) - c^in- l)(M-7i + l) , 



That the transmission is perfect can be shown by first re- 
alizing that in an interaction picture with respect to = 
Pi)/'^ ill which the diagonal elements of V vanish 
(this interaction picture will leave all entanglement properties 
unaffected as it is of direct sum form) we can replace 



Qit) 
Pit) 



cxp 



Vi 
-Vj 



(86) 



where Q and P are column vectors, by the complex notation 
Q -iP so that 



{Q~iP){t)=e'''''{Q^iP) 



(87) 



1 + 0.8086165(71 -l)-2/3 



(82) 



Clearly, for large n, the speeds approach a constant dependent 
on c. For the RWA interaction, the propagation velocity in- 
creases linearly with c. This is an attractive feature because 



Now we need to realize that Vj is a quantum mechanical rep- 
resentation of a rotation. This allows the evaluation of the 
matrix elements of e*^^*. In particular, we find that 



(e'^^*)iM = (sin(ct/2)) 



M-l 



(88) 
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so that one can generate an interchange between the first and 
the A/-th coordinate by waiting for a time t = tt/c (see Ref. 
10] for an analogous argument in spin chains). 

Without the assumption of the rotating wave approxima- 
tion, i.e., choosing the Hamiltonian i/spii„g to describe the 
time evolution the above simple argument fails and indeed it 
is not possible to tune the nearest neighbour couplings alone 
to generate perfect transfer of entanglement for AI > 2. How- 
ever, if one chooses the couplings as above and decreases the 
value of the constant c then, for a fixed distance, one can ob- 
tain arbitrarily good transfer efficiency at the expense of an 
increased delay time. This should not come as a surprise, as 
in the case of c — > the rotating wave approximation becomes 
exact as the terms that are neglected in the rotating wave ap- 
proximation are of order ©(c^). Therefore, for entanglement 
distribution over a fixed distance they will play a decreasing 
role as the time of arrival for the entanglement is of the order 
O(c-i). 

The case M = 2 is an exception, where one may realize an 
exact swap of the state of the 1-st to the 2-nd oscillator To 
show this, note that specific covariance matrix elements of the 
0-th and the 2-nd oscillator are given by 

7poP2(i) =-sinh(r)/i(t), 7go9„(i) = sinh(r)/i (t)(89) 
= sinh(r)/i(t), 7g2PoW = - sinh(r)gi (t). 

Considering the functions /i and gi as specified in eqs. i3Q\ . 
we find that there exists a real number c and a time t such that 
simultaneously 

= /iW=0, 5iW=0 (90) 



V = 



where Cij = c + Ac^.j is the position dependent coupling 
between the i-th and j-th oscillators, where Aaj is a real- 
ization of a random variable distributed according to a normal 
A^(0, Ac) distribution. For a chain of length 10, an average 
coupling constant of c = 0. 1 and an initial two mode squeezed 
state with squeezing parameter r = 0.8. Figure ^2 shows the 
ratio of the first maximum for the case of slightly perturbed 
couplings over the idealized case versus the perturbation size 
Ac. We observe that for Ac/c < 0.25, the achieved entangle- 
ment is greater than 95% of the degree of entanglement in the 
unperturbed case. Similar results apply for the RWA interac- 



can be satisfied. This is the case when we chose the real c 
such that there exist natural numbers k and I such that 

and t = Itt. Then, it follows that, as the 0-th oscillator is in- 
variant and the state of the 0-th and the 2-nd oscillator neces- 
sarily corresponds to a pure Gaussian state, the 1-st oscillator 
is necessarily decoupled from the other two. In this sense the 
state can be swapped from one oscillator to the other, while re- 
taining the entanglement with the 0-th oscillator Hence, one 
has a perfect channel for appropriate times. 



F. Sensitivity to random variations in the coupling 

In the preceding subsections we have discussed an ideal 
model in which all experimental parameters can be deter- 
mined perfectly. Any real experimental setup however will 
suffer small variations in parameters such as the coupling 
strength between neighbouring oscillators. In order to con- 
firm that the effects that have been found in this work can 
be observed in real experiments, we consider in the follow- 
ing the impact of random position dependent variations in the 
coupling strength between neighbouring oscillators. As an ex- 
ample, we consider an open chain with potential matrix 



(92) 



tion. These results indicate that sending quantum information 
along the chain is stable under perturbations. 

Similar considerations can be made for the spontaneous 
creation of entanglement which show that the results are much 
more sensitive to perturbations. Indeed, with the same specifi- 
cations as above we find the entanglement at the first maxima 
to be between a small fraction and twice the amount for the 
non-perturbed case when Ac/c = 0.5. This suggests that the 
experimental demands for the verification of the spontaneous 
creation of entanglement are considerably higher than for the 
distribution of entanglement. 



1 

1+Ci^2 -Cl,2 

-Ci,2 














^n—l.n 
'^n—l.n 1 H~ *^n— 1, 
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FIG. 1 1 : The ratio of the first maximum for the case of sHghtly per- 
turbed couplings over the idealized case versus the relative perturba- 
tion size Ac/c for a chain of length 30, r = 0.8 and c = 0.1. Each 
data point has been obtained as an average over 4000 realizations. 
The resulting curve is very well fitted by the function /(Ac/c) = 
1 + 0.02382Ac/c - 1.60481(Ac/c)^ + 1.59676(Ac/c)^ 



G. Other geometrical arrangements: Beamsplitters and 
interferometers 

So far we have studied only a linear chain of harmonic os- 
cillators through which quantum entanglement can be prop- 
agated. However, it might be interesting to consider more 
complicated structures which may be used as building blocks 
for more complicated networks, in principle, any arrangement 
coiTesponding to an arbitrary weighted graph. In this subsec- 
tion we will study briefly two possible extensions of the linear 
chain, namely a Y-shaped configuration which can be used 
for the generation of entanglement and a configuration resem- 
bling an interferometer. We show furthermore how such con- 
figurations may be switched on and off thereby controlling the 
transport of quantum information in such a structure. A more 

o 

o 

o 

o 

• o o o o o o 

o 

\ 

FIG. 12: This figure depicts a Y-shaped structure as it is described in 
the text. A single incoming arm consists of A/in oscillators (here 7) is 
connected to two outgoing arms each consisting of A/out oscillators 
(here 4). 



detailed discussion of such structures and their optimization 
will be presented elsewhere. The material in this subsection 
should merely serve as examples for possible alternative ways 
of creating and manipulating entanglement through propaga- 
tion in pre-fabricated structures. 

We begin by considering a chain in Y-shape which is shown 
in Figure One arm consisting of M\a oscillators is con- 
nected to two further arms each consisting of il/out oscillators. 
As usual we consider nearest neighbour interactions only and, 
for simplicity and the clearest demonstration of the effects, 
we restrict attention to the RWA interaction. We assume that 
the structure is initially in the ground state, i.e., at temperature 
T = 0. At time t = we perturb the first harmonic oscillator 
exciting it either to a thermal state characterized by covariance 
matrix elements 

T?!?! ~ Ipipi ~ ^ (93) 

for some z, or a pure squeezed state characterized by covari- 
ance matrix elements 

7gigi = l/7pipi = ^- (94) 

As an example we choose a coupling constant of c = 0.2 
and let the arms of the Y-shape contain A/out = 30 oscillators 
each while the base contains A/jn = 10 oscillators. In Figure 
1131 we present the results for the choice of a squeezed state 
with 7qigi = 10 = l/7pipi and a thermal state with ^q^q^ = 
10 = 7pipi- 

We observe that for an initial thermal state excitation no en- 
tanglement is ever found between the ends of the two arms of 
the Y-shape. This can be understood because a thermal state 
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'100 150 200 250 

Time 

FIG. 13: In a chain with Y-shape and nearest neighbour interaction 
of RWA type the first oscillator at the foot of the Y-shape is either 
excited to a squeezed state with ^qq — 10 = 1/ ^pp or a thermal state 
with ^qq = 10 = 7pp. The remaining oscillators are in the ground 
state. The perturbation propagates along the chain into both arms of 
the Y-shape. For an initial thermal state excitation no entanglement 
is ever found between the ends of the two arms of the Y-shape while 
entanglement is generated when the initial state is a squeezed state. 
The coupling constant is chosen as c = 0.2, the arms of the Y-shape 
contain 30 oscillators each while the base contains 10 oscillators. 
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is a mixture of coherent states, i.e., displaced vacuum states. 
If the system is initialized in the vacuum state it will evidently 
not lead to any entanglement in the RWA and therefore an 
initialization in a thermal state cannot yield entanglement ei- 
ther. On the other hand considerable entanglement is gener- 
ated when the initial state is a squeezed state. It is possible 
to optimize the generation of entanglement by adjusting the 
strength of the nearest neighbour couplings but this will be 
pursued elsewhere. These two observations are resembling 
closely optical beamsplitters which do not create entangle- 
ment from thermal state input but can generate entanglement 
from squeezed inputs (see Ref. |;35] for a comprehensive treat- 
ment of the entangling capacity of linear optical devices). 

Another interesting setup is shown in Figure [21 We will 
henceforth call this the interferometric setup. Let the num- 
ber of oscillators on the left (including the junction), up arm, 
down arm and on the right be AIl, Mjj, Mo and M^t respec- 
tively. If we prepare a two-mode squeezed state between a 
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FIG. 14: A diagram of the interferometric setup. 

decoupled oscillator and the leftmost oscillator of the inter- 
ferometric setup, then we are interested in how much entan- 
glement propagates through the setup depending on the prop- 
erties of the two arms. One may vary different parameters 
such as the length of one of the arms, the coupling strength or 
eigenfrequency of the oscillators in one arm. We will focus on 
how the change in eigenfrequency lo of the harmonic oscilla- 
tors in one of the arms affects the propagation of entanglement 
through the interferometric device. We change the eigenfre- 
quencies of the oscillators smoothly across one arm following 
cj,j = 1 + - 1) X min(i,A% + 1 - i)/{Mu/2) so that the 
oscillator half way through the arm has frequency uj. Figure 
[21 shows the logarithmic negativity between the decoupled 
oscillator and the last one in the interferometric configuration 
at the time t = 250 plotted against oj. The other parameters 
are A% = Md = 30, Ml = Mn. - 1 = 9 and c = 0.2. One 
clearly observes interference fringes in the frequency ui that 
are related to the effective path-length difference between the 
upper and the lower arm. The interference fringes do not have 
full amplitude and their amplitude is reduced for increasing 
u!. More sophisticated choices for the coupling parameters in 
the interferometric structure can improve on these imperfec- 
tions. This demonstrates that in interferometric structure the 
transmission through the device will be strongly influenced by 
changes of the properties of one arm of the structure. 

This shows that more complicated structures such as the 
Y-shape or the interferometric may be used to create entan- 





FIG. 15: The logarithmic negativity N{uj) between the decoupled 
oscillator and the last one in the interferometric configuration at the 
time t — 250 plotted against u). We have chosen Ml = Mr — 1 = 9, 
Mu = Md = 30 and c = 0.2 in the RWA interaction. One clearly 
observes interference fringes. 



glement from an initially unentangled system and transport it. 
There is a distinct analogy here to quantum optical networks 
which might be used for information processing either em- 
ploying the polarization degree of freedom or as we did here 
the excitation number degree of freedom. This suggests that 
one could construct similar 'hardwired' networks on the level 
of interacting quantum systems that could then perform cer- 
tain quantum information processing or communication tasks. 
This might involve structures such as the Y-shape presented 
here but may also implement structures such as interferome- 
ter structures shown in Figure^Jor multi-input devices. 

If one were to consider hardwired structures, then it would 
be necessary to devise methods by which these structures 
could be switched on and off. Here we explore two possibili- 
ties. Firstly, one might change the coupling strength cjunction 
of the oscillator at the three way junction in the Y-shape. 
Apart from the obvious fact that they remain disentangled for 
Cjunction = (i.e., uncouplcd), we find the first maximum for 
the entanglement decreases to roughly half the value for large 
coupling strength Cjunction = 0.8. This is shown in FigureFTSl 
A further increase of the coupling strength to cjunction = 5 
does not lead to further significant change. A different ap- 
proach would be to change the eigenfrequency or the mass of 
the junction oscillator while keeping the coupling strength the 
same as the other oscillators. Indeed, if we increase the eigen- 
frequency (^Junction, the entanglement can be reduced to an ar- 
bitrarily small amount for both RWA and spring-interactions. 
A decrease of the eigenfrequency is less efficient but would 
also allow a significant reduction of the amount of entangle- 
ment generated in the device. FigurefTTldemonstrates these ef- 
fects achieved by changing the eigenfrequency of the junction 
oscillator in the Y-shape. It should be noted that the depen- 
dence of the logarithmic negativity with the eigenfrequency 
exjunction IS almost perfectly fitted by a Lorentzian line shape. 
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FIG. 16: The degree of entanglement in the Y-structure in terms 
of the coupling strength cjunction of the oscillator at the three-way 
junction. The parameters are chosen identically as in Figure fTsl 




FIG. 17: By increasing and decreasing exjunction respectively in the 
Y-shape we observe a noticeable change in the amount of entangle- 
ment that is generated in the device. The other parameters are chosen 
as in Figure [TbI 



The above examples suggest, that it is possible to switch 
on and off pre-fabricated devices such as the Y-shape shown 
above by adjusting either the coupling strength (decreasing it, 
i.e., approach decoupling) or the eigenfrequency (increasing 
it). Such a manipulation of the junction oscillator dictates the 
quantum information flowing through the junction. An exam- 
ple for a possible implementation of such a switch from an op- 
tical setting would be coupled optical cavities which are filled 
with atoms. Laser irradiation of these atoms would then lead 
to a shift of the resonance frequency of the cavity which cor- 
responds to a change in the eigenfrequency of a harmonic os- 
cillator in the above examples. In this way, individual cavities 
might be decoupled. A detailed study of such a scheme will be 
presented elsewhere. Finally, we would like to briefly mention 
an analogy with the monotonic and non-monotonic behaviour 



of the efficiency. We find that by increasing the mass of the 
junction oscillator in the beamsplitter configuration, we obtain 
monotonically decreasing entanglement between the two ends 
for the RWA interaction whereas the Hooke's law interaction 
produced non-monotonic behaviour. 



VI. CONCLUSIONS AND DISCUSSIONS 

We have investigated the entanglement dynamics of sys- 
tems of harmonic oscillators both analytically and numeri- 
cally. Particular attention has been paid to harmonic oscil- 
lators coupled by springs (Spring) and to harmonic oscilla- 
tors with a linear coupling in a rotating wave approximation 
(RWA) as it is appropriate in a quantum optical setting. After 
an introduction to the mathematical formalism and the deriva- 
tion of the analytical solutions for the equations of motion for 
these interactions we then investigated several possible sce- 
narios. We considered the generation of entanglement with- 
out detailed local control of individual systems. This was 
achieved by first switching off any interaction between the 
oscillators, cooling them to near the ground state and subse- 
quently switching on the coupling suddenly. Surprisingly, en- 
tanglement will be generated over very large distances which 
is in stark contrast to the entanglement properties of the sta- 
tionary ground state of a harmonic chain where only nearest 
neighbours exhibit entanglement fl^ . We have also demon- 
strated that a linear chain of harmonic oscillators is able to 
transport quantum information and quantum entanglement for 
various types of nearest neighbour coupling. For position 
independent nearest neighbour coupling we observe that the 
transmission efficiency is a non-monotonic function in the 
coupling strength for Hooke's law coupling while it is mono- 
tonically decreasing for the RWA coupling. In both cases this 
suggests that it is advantageous to transmit entanglement in 
smaller portions rather than large units. But due to the rapid 
decline in efficiency with the spring interaction for very small 
r, one should avoid sending in too smaller r. The propaga- 
tion speed for the quantum entanglement has been provided 
analytically. For the above effects we have studied their sen- 
sitivity to random variations in the coupling between the os- 
cillators and to finite temperatures. 

Finally we have proposed more complicated geometrical 
structures such as Y-shapes and interferometric setups that al- 
low for the generation of entanglement in pre-fabricated struc- 
tures without the need for changing any coupling constants. 
We have also shown that these structures may be switched 
on and off by changing the coupling of only a single har- 
monic oscillator with its neighbours. This suggests the pos- 
sibility for the creation of pre-fabricated structures that may 
be 'programmed' by external actions. Therefore quantum in- 
formation would be manipulated through its propagation in 
these pre-fabricated structures somewhat analogous to modern 
micro-chips and as opposed the most presently suggested im- 
plementations of quantum information processing where sta- 
tionary quantum bits are manipulated by a sequence of exter- 
nal interventions such as laser pulses. 

All these investigations were deliberately left at a device in- 
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dependent level. It should nevertheless be noted that there are 
many possible realizations of the above phenomena. These 
include nano-mechanical oscillators fllll . arrays of coupled 
atom-cavity systems, photonic crystals, and many other re- 
alizations of weakly coupled harmonic systems, potentially 
even vibrational modes of molecules in molecular quantum 
computing ll36ll . A forthcoming publication will discuss de- 
vice specific issues of such realization as well as well as im- 
proved structures (including novel topological structures as 
well as changes of their internal structure such as di-atomic 
chains) that allow for better performances with less experi- 
mental resources. We hope that these ideas may lead to the 
development of novel ways for the implementation of quan- 
tum information processing in which the quantum information 



is manipulated by flowing through pre-fabricated circuits that 
can be manipulated from outside. 

We acknowledge discussion with D. Angelakis, S. Bose, 
D. Browne, M. Christandl, M. Cramer, J. DreiBig, J.J. HalU- 
well, T. Rudolph, M. Santos, V. Vedral, S. Virmani, and R. 
de Vivie-Riedle. This work was supported by the Engineer- 
ing and Physical Sciences Research Council of the UK, the 
ESF programme "Quantum Information Theory and Quan- 
tum Computing", the EU Thematic Network QUPRODIS and 
the Network of Excellence QUIPROCONE, the US Army 
(DAAD 19-02-0161), the EPSRC QIP-IRC, and the Deutsche 
Forschungsgemeinschaft DFG. MBP is supported by a Royal 
Society Leverhulme Trust Senior Research Fellowship. 



[1] J.I. Cirac, P. Zoller, H.J. Kimble, and H. Mabuchi, Phys. Rev. 

Lett. 78, 3221 (1997). 
[2] L.-M. Duan and H.J. Kimble, Phys. Rev. Lett. 90, 253601 

(2003). 

[3] X.L. Feng, Z.M. Zhang, X.D. Li, S.Q. Gong, Z.Z. Xu, Phys. 

Rev. Lett. 90, 217902 (2003). 
[4] D.E. Browne, M.B. Plenio and S.F. Huelga, Phys. Rev. Lett. 91 

067901 (2003). 

[5] N. Khaneja and S.J. Glaser, Phys. Rev. A 66, 060301 (2002). 

[6] S. Bose, Phys. Rev. Lett. 91, 207901 (2003). 

[7] V. Subrahmanyan, LANL e-print quant-ph/0307135 

[8] M. Chi'istandl, N. Datta, A. Ekert, A.J. Landahl, LANL e-print 

quant-ph/0309131 

[9] T. Osborne and N. Linden, LANL e-print quant-ph/0312141 
[10] J. Eisert and M.B. Plenio, Int. J. Quant. Inf. 1, 479 (2003), 

LANL e-print quant-ph/03 1 207 1 
[11] J. Eisert, M.B. Ple nio, S. Bose and J. Hartley, LANL e-print 

|quant-ph/03111I3| 
[12] K. Audenaert, J. Eisert, M.B. Plenio, and R.F. Werner, Phys. 

Rev. A 66, 042327 (2002). 
[13] T.A. Brun and J.B. Hartle, Phys. Rev. D 60, 123503 (1999). 
[14] J.J. Halhwell, Phys. Rev. D 68, 025018 (2003). 
[15] M.L. Roukes, Phys. World 14, 25 (2001). 
[16] K. Schwab, E.A. Henriksen, J.M. Worlock, and M.L. Roukes, 

Nature 404, 974 (2000). 
[17] M.M. Wolf, F. Verstraete, and J.I. Cirac, Phys. Rev. Lett. 92, 
087903 (2004). 

[18] M. Hein, J. Eisert, and H.J. Briegel, LANL e-print 

quant-ph/0307130 
[19] D. Schlingemann, Quant. Inf. Comp. 2, 307 (2002). 
[20] R. Raussendorf, D.E. Browne, and H.J. Briegel, Phys. Rev. A 

68, 022312 (2003). 



[21] R. Simon, E.C.G. Sudarshan, and N. Mukunda, Phys. Rev. A 
36, 3868 (1987). 

[22] K. Audenaert, M.B. Plenio, and J. Eisert, Phys. Rev. Lett. 90, 
027901 (2003). 

[23] K. Zyczkowski, P. Horodecki, A. Sanpera, and M. Lewenstein, 

Phys. Rev. A 58, 883 (1998). 
[24] G. Vidal and R.F Werner, Phys. Rev. A 65, 032314 (2002). 
[25] J. Eisert and M.B. Plenio, J. Mod. Opt. 46, 145 (1999). 
[26] J. Eisert, PhD thesis, Potsdam, February 2001. 
[27] D. Guilini, E. Joos, C. Kiefer, J. Kupsch, I.O. Stamatescu, and 

H.D. Zeh, Decoherence and the Appearance of a Classical 

World in Quantum Theory, (Springer, Berlin-Heidelberg-New 

York, 1996). 

[28] W.H. Zurek, Rev. Mod. Phys. 75, 715 (2003). 

[29] H.P. Breuer and F. Petruccione, The Theory of Open Quantum 
Systems (Oxford University Press, Oxford, 2002). 

[30] U. Zurcher and P Talkner, Phys. Rev. A 42, 3267 (1990). 

[31] Royal Society Mathematical Tables Volume 7, Bessel Functions 
Part 3, Zeros and Associated Values, F.W.J. Olver ed. (Cam- 
bridge University Press, Cambridge, 1960). 

[32] Note that a similar non-monotonic behaviour of efficiencies can 
be observed in different situations such as in the efficiency of 
transferring entanglement from squeezed light to two qubits 
I33I1 or the generation of entanglement from thermal light I34I1 . 

[33] B. Kraus and J.I. Cirac, Phys. Rev. Lett. 92, 013602 (2004). 

[34] M.B. Plenio and S.F Huelga, Phys. Rev. Lett. 88, 197901 
(2002). 

[35] M.M. Wolf, J. Eisert, and M.B. Plenio, Phys. Rev. Lett. 90, 
047904 (2003). 

[36] CM. Tesch and R. de Vivie-Riedle, Phys. Rev. Lett. 89, 157901 
(2002). 



